
function [p q] = implicit_theta(dt,beta)
global U  h  g  p  q  nt dx x h1 p_global dx1 dx2 xrange;;

dt2=1/dt;
A = zeros(x,x);

A(1,2) = U*0.5*dx1(1);
A(1,x) = -U*0.5*dx1(1);
% reference: http://onlinelibrary.wiley.com/doi/10.1111/j.2153-3490.1970.tb01933.x/pdf
for k=2:x-1
    r = dx(k)/dx(k-1);
    A(k,k) = -U*(1-r*r)/(dx(k)*(1+r));
    A(k,k+1) = U/(dx(k)*(1+r));
    A(k,k-1)= -U*r*r/((dx(k)*(1+r)));
    
%     A(k,k-1) = -U/(dx(k-1)+dx(k));
%     A(k,k+1) =  U/(dx(k-1)+dx(k));
end

A(x,x-1) = -U*0.5*dx1(x-1);
A(x,1)   =  U*0.5*dx1(x-1);

B = zeros(x,x);

B(1,1) = -2*h*dx2(1);
B(1,2)= h*dx2(1);
B(1,x)=h*dx2(1);
    
for k=2:x-1;
    r = dx(k)/dx(k-1);
    dm = dx(k)*dx(k-1)*(1+r);
    
    B(k,k+1) = 2*h/dm;
    B(k,k-1) = 2*h*r/dm;
    B(k,k) = -2*h*(1+r)/dm;
    
    
%     B(k,k+1)= h*2*(1/(dx(k)*(dx(k-1)+dx(k))));
%     B(k,k-1)= h*2*(1/(dx(k-1)*(dx(k-1)+dx(k))));
%     B(k,k) = -2*h/(dx(k-1)*dx(k));
end

B(x,x) = -2*h*dx2(x-1);
B(x,x-1)= h*dx2(x-1);
B(x,1)=h*dx2(x-1);
L=zeros(2*x,2*x);
L(x+1:2*x,x+1:2*x)=A;
L(1:x,1:x)=A;
L(1:x,x+1:2*x)=B;

for k=1:x
    L(x+k,k)=g;
end


I = eye(2*x);
sol(1:x,1) =p;
sol(x+1:2*x,1) =q;
%p_global(:,1)=p;

for n=2:nt+1;   
    rhs = (dt2*I - (1-beta)*L)*sol;
    sol=(dt2*I+beta*L)\rhs;    
    p = sol(1:x,1);
    q= sol(x+1:2*x,1);
%p_global(:,n)=p;
    if rem(n,100)==0
        
        refreshdata(h1,'caller') % Evaluate p in the function workspace
        drawnow
    end
    
end
%display('Completed Successfully');